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Abstract. This paper deals with the time evolution in the 
matter era of perturbations in Friedman-Lemaitre models 
with arbitrary density parameter £1, with either a zero 
cosmological constant, A = 0, or with a non-zero cosmo- 
logical constant in a spatially flat Universe. Unlike the 
£^ classical Eulerian approach where the density contrast is 
^ expanded in a perturbative series, this analysis relies in- 
— 5 stead on a perturbative expansion of particles trajectories 
in Lagrangian coordinates. This brings a number of ad- 
vantages over the classical analysis. In particular, it en- 
C*~) ables the description of stronger density contrasts. Indeed 
' _| the linear term is the famous Zel'dovich approximate so- 
lution (1970). This approach was initiated by Moutarde et 
al. (1991), generalized by Bouchet et al. (1992), and fur- 
ther developed by many others. We present here a system- 
atic and detailed account of this approach. We give ana- 
H lytical results (or fits to numerical results) up to the third 
O jnrder (which is necessary to compute, for instance, the four 
q point spatial correlation function or the corrections to the 
5^ linear evolution of the two-point correlation function, as 
c/2 well as the secondary temperature anisotropies of the Cos- 
& mic Microwave Background). We then proceed to explore 
the link between the lagrangian description and statistical 
measures. We show in particular that Lagrangian pertur- 
bation theory provides a natural framework to compute 
the effect of redshift distortions, using the skewness of the 
density distribution function as an example. Finally, we 
show how well the second order theory does as compared 
to other approximations in the case of spherically sym- 
metric perturbations. We also compare this second order 
approximation and Zel'dovich solution to N-body simula- 
tions in the description of large-scale structure formation 
starting from a power law (n = —2) power spectrum of 
Gaussian perturbation. We find that second order theory 
is both simple and powerful. 
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1. Introduction 

The observed large scale structures revealed by galaxy cat- 
alogs are usually supposed to have arisen due to gravita- 
tional instability acting on small initial perturbations in 
the context of an expanding Universe. The corresponding 
dynamics may in principle be described either from an 
Eulerian or Lagrangian point of view. With the notable 
exception of Zel'dovich (1970) approximate solution (and 
numerical simulations), analyses have up to very recently 
been largely dominated by the use of the Eulerian ap- 
proach. Actually, most Eulerian studies concentrated on 
the perturbative regime, when density contrasts and ve- 
locities are small, and calculations are tractable. By prin- 
ciple, thus, such an approach is limited to the description 
of weak density contrasts. This is not the case of the La- 
grangian approach, as is exemplified by the case of one- 
dimentional perturbations when Zel'dovich approximation 
is in fact exact up to the singularity, when the density con- 
trast becomes infinite! 

In order to contrast the two points of view, we 
start by recalling the classical Eulerian perturbative ap- 
proach. Then we review some recent developments con- 
cerning the Lagrangian perturbative approach introduced 
by Moutarde et al. (1991). In this paper, as in Bouchet 
et al. (1992), we study directly the Newtonian perturba- 
tive approach, and extend this work in several ways. In the 
next section (§2), we establish notations and give some de- 
tails of the derivation which was, due to lack of space, only 
sketched in the letter by Bouchet et al. (1992). In keeping 
with Zeldovich's spirit, we focus on the fastest growing 
modes of irrotational flows. We also provide quite accu- 
rate approximations to the third order solutions (which 
cannot be expressed exactly, for SI ^ 1, by using simple 
mathematical functions). Similarly we give expressions up 
to third order in the flat A / case. In §3 we relate 
statistical indicators to the properties of the initial con- 
ditions. This is applied to the real space-redshift space 
mapping of low order statistical properties. Section 4 is 
devoted to a systematic comparison of the Lagrangian 
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second order solution with other approximations, in the 
spherically symmetric case whose exact evolution is well- 
known. First and second order Lagrangian approximations 
are then compared to N-body simulations in the case of 
Gaussian initial conditions with a power-law power spec- 
trum (n = —2) in §5. We conclude in §6 by a discussion. 

1.1. Classical Eulerian Approach 

Our aim here is to recall the basic steps of the pertur- 
bative Eulerian approach. The same conceptual steps are 
involved in the Lagrangian case. But the simple change of 
"small parameter" governing each of the two expansions 
(i.e. the density contrast for the Eulerian approach ver- 
sus the displacement field for the lagrangian one) has far 
reaching consequences, in terms of the range of validity of 
each approach, and more generally in terms of the class of 
problems tractable analytically. A more detailed account 
of the Eulerian theory can be found for instance in Peebles 
(1980, hereafter LSS). 

Consider perturbations which may be described in the 
Newtonian approximation. It is convenient to use comov- 
ing coordinates, 

2 dx 

r = ax, and n = ma ——, 
1 dt ' 

where m and a stand respectively for the particles mass 
and the scale factor of the metrics of the Friedman- 
Lemaitre background model. The motion and field (Pois- 
son) equations then read 



dp 
~dt 



-f- = — V0, V 2 <j> = 4irGa 2 p8 



(1) 



if 6 = p/~p — 1, and ~p is the mean mass density. 

The kinetic theory equation governing the evolution of 
the system then obtains by requiring mass to be conserved. 
If f(x,p,t) stands for the probability density of finding 
at time t a collisionless particle within the infinitesimal 
volume dxdp, Liouville theorem yields the Vlasov (or col- 
lisionless Boltzmann) equation 



dtf- 



P 



V/-mVf 9 p / = 0, 



where partial derivatives over t are denoted by dt- Taking 
velocity moments of this mean-field equation then leads to 
an infinite hierarchy. In particular, the macroscopic den- 
sity, p = ma~ 3 J dp f = ~p(l + 8), and macroscopic ve- 
locity, v = f dpp f/(ma 2 f dpf), satisfy the conservation 
equation 



apd t 8 + V(pv) = 
and 

o t 20 + 2 o = - 

a a 



(2) 



V-[(1 + <S)V<? 



1 



rd a dp [(1 + 8) (p 

aPf))\ 



(3) 



In the pressureless case, ((p a — mav a )(p/3 — mavp)) = 0; 
thus (p a pfj) l(m 2 a 2 ) = v a vp, which closes the hierarchy. 
Alternatively, one could start from the perfect fluid equa- 
tion (LSS §5, p.47). 

Perturbative solutions are then obtained by means of 
an iterative procedure. Let us write the density constrast 
as 

S = £ 8W+s 2 8^+e 3 8^ + 0(e 4 ), 

where e is just a book-keeping device, and j < 1. First, 
keeping only the terms of order e, one obtains 

5t ^(i) + 2^^ 1 ) = 47rGa 2 ^ 1 ), 
a 

which governs the linear evolution. 

In order to solve this equation, one now has to specify 
the time variation of the scale factor. In the absence of 
a cosmological constant, and in the matter era when ~p oc 
a -3 , Friedman equation can be rewritten as follows 



d t a 
a 



8jrGp 



8ttGp 



a 2 R 2 



l + (n -i_l)_ 

do. 



with the subscript corresponding for instance to today. 
In the Einstein-De Sitter case, when SI = 1, the solution 
is a oc t 2 l 3 . The general solution is then a linear superpo- 
sition of growing and decaying modes, D a and Dj, 



8 (1) = k a L 
with 

D -t 2/s 



kh Dh 



id D h 



(4) 



The corresponding velocity field (actually its divergence) 
is deduced from the conservation equation (2), and the 
constants k a and ki, are determined by the initial condi- 
tions for the density (or the gravitational potential) and 
the velocity fields. 

For an open model with SI < 1, the parametric solution 
of Friedman equation is 



a = j4(cosh r) — 1), 
with 

4ttG 



id t = 5(sinhry — 1), 



A 



-pa d \R\* 



a 



3 r 2 11-Slo 1 ! 



B = A\R\ 



(5) 



It is then convenient to use y 
able, in which case 



SI 1 — 1 for time vari- 



D a 



(i + y) 112 



3/2 



(6) 
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Finally, the closed model case obtains by transforming the 
equations (rj, A, B, R, y) — ► (irj, —A, iB, iR, —y). 

If one now keeps also the 0(s 2 ) terms in eq. (3), one 
gets 

dt*6 + = 4irGa 2 p8(l + 8) + V8-^ + \d a dp(v a v p 

a a 1 a 1 

The solution in the Einstein-De Sitter case is found by 
using the first order solution (4) 
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47r" 



8 (1 JA_. 



1 



567T 2 
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where A stands for —<f)/(G~pa 2 ). As was noted by Martel 
and Freudling (1991), the equation (3) does not appear 
separable for SI ^ 1. But they found numerically a so- 
lution by assuming SI 1 ' 2 k, £1. We shall in the following 
derive the exact analytical solution found by Bouchet et 
al. (1992) for arbitrary values of SI, in the absence of any 
cosmological constant. We shall also give an approxima- 
tion of the third order solution, including the A ^ case, 
in terms of usual mathematical functions. 

1.2. Developments of the Lagrangian Approach 

Zeldovich (1970) proposed to overcome some of the dif- 
ficulties of the (pressureless) Eulerian approach by using 
the linear theory in terms of Lagrangian coordinates. The 
primary object of the analysis is then the trajectory of a 
"particle" instead of the density contrast. A fluid element 
(a particle) is indexed by its unperturbed Lagrangian co- 
ordinate q. Its comoving Eulerian position at time t, x(i), 
is connected to q by a displacement field 



x = q + &(t,q). 



(7) 



Zeldovich kinematic approximation amounts to assume 
that the displacement field is simply proportional to the 
initial displacement field, \P l (q), 

*(t,q) = g 1 (t)* i (q), 

which is indeed a self-consistent irrotational solution of 
the linearized equations of motion, provided gi(t) is given 
by the linear growth rate (4). In effect, this is a ballistic 
approximation where gravitational acceleration is ignored; 
the movement is simply inertial. 

The idea was to use this ansatz, even when the density 
contrast 6 becomes large. It lead to the development of 
pancake theory. Indeed, the density constrast is given by 
the determinant, J, of the Jacobian of the transformation 
from x to q, 



1_ 
7 



1 



with J = | detX>|, V being the tensor of deformation, 



V 



dx 
dq_ 



2 + T, T 



dq _ 



(8) 



where X is the identity matrix. If we call A; the eigenvalues 
of T/}oo , then J = J^[? =1 (l + g\ (t)A{), and there is a pan- 
cake forming perpendicular to the principal axis ofD with 
the largest negative eigenvalue (if any). The spatial struc- 
ture of perturbations near maxima of the density field (and 
)the connection to angular momentum of galaxies and clus- 
ters) was thoroughly investigated by Doroshkevich (1970). 
Further developments along these lines are reviewed by 
Zeldovich (1978) and Doroshkevich, Shandarin, and Saar 
(1978). 

Of course, as stated by Zeldovich, "the analytic evalua- 
tion of the error is extremely difficult" . This lead to many 
comparisons with the exact dynamics, either with simu- 
lations (e.g., Doroshkevich et al. 1980) or rigorous Eule- 
rian perturbative theory (Grinstein and Wise 1986). The 
approximation turns out to be amazingly good, at least 
when the initial field is smooth enough. Indeed, it is widely 
used today, e.g., to predict the weakly non-linear evolu- 
tion of the moments of the one point probability distri- 
bution function of the density field (Betancort-Rijo 1991; 
see also Hoffman 1987 for the variance only), or of the 
distribution itself (Kofman et al. 1993, Padmanabhan & 
Subramanian 1993), or in reconstruction methods to re- 
cover the "initial conditions" from present day observa- 
tions (e.g., Nusser et al. 1991, Nusser and Dekel 1992, 
Gramman 1992, Lachieze-Rey 1993a). In fact, it is even 
used rather often to address questions where it is not ap- 
propriate. 

The success of Zeldovich approximation brought about 
many attempts to do better by correcting its shortcom- 
ings. A limitation of Zeldovich approximation is that, be- 
ing linear (i.e., neglecting the gravitational acceleration 
induced by the changing density in the course of evolu- 
tion), pancakes, once formed, indefinitely thicken in the 
absence of a restoring force (and also collapse times are 
overestimated). This can be vividly illustrated by examin- 
ing Figure 1 which shows particle trajectories in a simple 
toy model. Following Moutarde et al. (1991) and Buchert 
(1993) this model is the superposition of 2 orthogonal sine 
waves with the same amplitude and wavelength, in a flat 
expanding Universe. This toy model is implemented in a 
square box, the axes of which are parallel to the waves 
and their comoving length, equal to the wavelength, is set 
to unity. 64 2 particles are initially laid down on a regu- 
lar lattice and at a = 1 their positions and velocities are 
perturbed, with an initial displacement field, IP* given by 




t sin(27ra;), 



(9) 



£sin(27rt/), 6i = ±, 



and (da;/da)i = \P l . The numerical simulation of the exact 
evolution was done with a PM code using a 128 2 grid. 
With such initial conditions, particles move toward the 
center of the box, creating an overdensity and a potential 
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Fig. 1. Trajectories of particles for 2 in-phase, perpendicular, sine waves. The exact trajectories computed by direct numerical 
simulation are shown at the top left. The top right panel shoes the trajectories in Zeldovich approximation. The bottom panel 
shows the trajectories obtained by using the frozen flow approximation (left), and the second-order Lagrangian approximation 
(right). 



well in which they will oscillate. The figure shows x(a) 
for the particles initially closest to y = 1/2. According 
to Zeldovich approximation a particle with coordinate q 
before the perturbation has an Eulerian coordinate 



x(a,q) 



(10) 



at expansion factor a. This results in straight line tra- 
jectories, with an overestimate of the collapse time, and 



artificial post-collapse thickening of the non-linear struc- 
tures. 

A popular remedy to the second shortcoming of Zel- 
dovich approximation is the "adhesion" approximation in- 
troduced by Gurbatov and Sai'chev (1981, 1984) which 
forces particles orbits not to cross by introducing an artifi- 
cial viscosity. This leads to Burger's equation of non-linear 
diffusion (Burgers 1940, 1974). It was later developed by 
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many others, in particular Gurbatov, Sai'chev, and Shan- 
darin (1985, 1989), and Shandarin (1987); see the review 
by Shandarin and Zeldovich (1989) for further references. 
It has mainly been used to study the very large scale struc- 
tures: for instance, Nusser and Dekel (1990) used this ap- 
proach to study numerically their filamentary aspect in 
scale-invariant models, while Weinberg and Gunn (1990) 
focused on CDM initial conditions. For recent and more 
complete references, see Sahni et al. 1994 and Vergassola 
et al. 1994). 

More recently, Matarese et al. (1992) suggested to im- 
prove Zeldovich approximation by replacing it with a 
"frozen flow" approximation. In Zeldovich solution, par- 
ticles always keep the same velocities; these velocities are 
in effect "frozen" . Matarese and its colleagues proposed 
instead to freeze the initial velocity field: particles move 
along streamlines computed (once and for all) from the 
initial potential; consequently, they never even reach the 
singularities. The figure 1 shows the corresponding tra- 
jectories in our toy model. Clearly, the pre-shell-crossing 
trajectories are not well described, but the approximation 
does capture some features of the asymptotically late evo- 
lution. A further improvement along this line of thought 
has just been proposed by Bagla and Padmanabhan (1994) 
who propose to freeze instead the initial potential (which 
is linearly conserved). 

A much more ambitious program would be to try 
to find exact solutions to the full Lagrangian equations 
(Buchert and Gotz 1987), maybe restricted to specific ini- 
tial conditions. The only tractable case so far (Buchert 
1989) is when the motion is locally one-dimensionnal, 
when two eigenvalues are (locally) zero. This is not too 
surprising, since Zeldovich approximate solution is well- 
known to be in fact exact for a purely one-dimensionnal 
collapse. More recently, Lachieze-Rey (1993) found a new 
class of formal solutions. 

To our knowledge, the apparently obvious idea of us- 
ing a rigorous Lagrangian perturbative expansion is due 
to Moutarde et al. (1991). These authors, in the course of 
a systematic study of smooth and isolated perturbations, 
first found numerically that the evolution of three sine 
waves leads, at shell crossing (when orbits intersect), to a 
scale-invariant density profile; their results were thus ex- 
tending the two-dimensionnal result of Alimiet al. (1990). 
They reasoned that this could be analytically confirmed 
by a Lagrangian approach, in the same spirit that lead Zel- 
dovich to his approximation, since this power-law density 
profile builds up before shell crossing. They performed a 
Lagrangian perturbative expansion up to the third order, 
in the SI = 1 case. Their determination of the spatial part 
of the solution, and of the growth rates up to the third or- 
der, was only performed explicitly for the initial conditions 
of their simulations. They found exact agreement with the 
numerical density profile near the collapse time up to a 
density contrast of order 50! Figure 1 shows the trajecto- 
ries for our toy model with sine waves as computed with 



the second-order lagrangian approximation (hereafter re- 
ferred to as L2). Indeed, as long as particle trajectories 
do not cross, even at very late stages, their trajectories 
are correctly described by the second order approxima- 
tion. Consider for instance the particles initially close to 
x = 0, their position up to a « 80 are well predicted. At 
the same time, trajectories of particles initially closer to 
x = 1/2 remain correct up to the time of collapse and 
therefore the time of first shell crossing is more correctly 
estimated by L2 than by Zeldovich approximation. On 
the other hand, Zeldovich approximation gives better re- 
sults after shell crossing because particles go apart from 
the point of collapse more slowly than in L2. One can see 
that at a = 80 structures are completely washed out in 
the L2 approximation. However, one can estimate that, 
as far as a ps 30, a bit after shell crossing, L2 is closer to 
the simulation than Zeldovich. Thus figure 1 clearly shows 
what to expect for the various approximations: the frozen 
flow captures some features of the asymptotically late be- 
havior, while L2 improves the description of the weakly 
non-linear phase. 

The work of Moutarde et al. (1991) was soon general- 
ized in two ways. On one hand, Buchert (1992) explored 
the perturbative solutions for arbitrary initial conditions 
(in the SI = 1 case). He found for instance that Zeldovich 
solution is a sub-class of the general first order solutions, 
when peculiar velocities and accelerations are required to 
be parallel (which is true for irrotational flows at late 
times, when the growing mode dominates, as was explic- 
itly assumed by Zeldovich from start). He also explored 
the first order solutions for rotational perturbations. 

On the other hand, Bouchet et al. (1992) explored a 
more restricted class (the same class of solution as Zel- 
dovich's, namely the fastest growing part of an irrota- 
tional flow), but gave the second-order exact solution for 
Friedman-Lemaitre models of arbitrary density parameter 
SI, provided there is no cosmological constant (the first or- 
der, i.e., Zeldovich solution for this case has been known 
for long, as in the Eulerian case, see Shandarin 1980; one 
even knows the behavior with a cosmological constant, 
Bildhauer et al. 1992). And as an illustration of the power 
of the approach, Bouchet et al. (1992) made a first connec- 
tion with measurable statistical quantities, by computing 
the value of the skewness of the density field probability 
distribution function (PDF). For SI = 1, it is a constant, 
S3 = 34/7, times the square of the variance of the PDF, 
as was already found by Peebles (1980). They showed that 
S3 depends extremely weakly on SI (see below). And that 
S3 is barely affected by a real space-redshift space map- 
ping. This is brought further by Juszkiewicz et al. (1993; 
see also Juszkiewicz and Bouchet 1991) who consider the 
effect of smoothing, and by a companion paper by Hivon 
et al. (1994) who consider the effect of smoothing in red- 
shift space. A recent summary of most of these results may 
be found in Bouchet and Juszkiewicz (1993). 
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Quite recently, Buchert and Ehlers (1993) and Buchert 
(1993) explore (respectively) the generic second and third 
order solutions in the £1 = 1, A = case, while Lachieze- 
Rey (1993) propose a tensorial reformulation of the equa- 
tions of Bouchet et al. (1992) which naturally includes the 
curl-free and divergence-free part of the vector field. In- 
terestingly, Lachieze-Rey found a special tensorial solution 
which is purely local. Unfortunately, it is not a general so- 
lution. Indeed, an alternative view of the problem is given 
by its relativistic generalization which was reviewed by 
Ellis (1971). It was then found that a quantity called the 
magnetic part of the Weyl tensor vanishes in the Newto- 
nian limit (when keeping the highest order term in a l/c 
expansion). The equations are then quite simple and, as- 
suming this holds at higher orders, Barnes and Rowlingson 
(1989), Matarese, Pantano, and Saez (1993), Bertschinger 
(1993), Croudace et al. (1994) obtained closed sets of local 
Lagrangian perturbative equations. 

This generated a lot of excitement, as it seemed to 
magically evade the non-local character of the perturba- 
tive Newtonian approach. Bertschinger and Jain (1994) 
proceeded to deduce that, contrary to Zeldovich pancak- 
ing picture, typical perturbations would rather collapse as 
filaments. Quite recently, though, it has been realized that 
this is incorrect, since at the post-newtonian level (i.e. at 
higher orders of a l/c expansion) the magnetic part of the 
Weyl tensor is non-vanishing and should be included in 
a non-linear treatment of the post-Newtonian limit (Kof- 
man and Pogosyan 1994; see also Bertschinger and Hamil- 
ton 1994). Indeed, a correct l/c expansion of the equations 
of General Relativity does lead to the same equations as 
in the standard perturbative Newtonian theory, implying 
that the Lagrangian evolution is not purely local. 

2. Lagrangian Perturbation Theory 

The Lagrangian perturbative approach introduced by 
Moutarde et al. (1991) proceeds in quite a parallel fashion 
to the Eulerian case recalled in the introduction. As be- 
fore, we use the Newtonian approximation with comoving 
coordinates, but following Doroshkevich et al. (1973) we 
replace the standard cosmological time by a new time r 
defined by 

dT(xa~ 2 dt. (11) 

In that case the motion and field equations now read 

i = -V*$, A x $ = /3(t)6, (12) 

with no x term. In this equation, dots denote derivatives 
with respect to r, and /?(t) = AitGa(J>a 3 ). By choos- 
ing the proportionality constant in the definition (11) to 
be — l,j4/|i?|, or — when, respectively, £1 is equal, 

smaller or greater than 1, [3 is simply given by 



r 2 + fc(Sl)' 



with 

k(Q = 1) = 

k(n < i) = -i 
k(n > i) = +i. 

While t = t' 1 ! 3 , when SI = 1, one has r = |1 - Sll" 1 / 2 
otherwise. 

Now we wish to follow the particle trajectories instead 
of the density contrast, by using the mapping (7). The 
Jacobian of the transformation from x to q, permits to 
express the requirement of mass conservation simply as 

p(x)Jd 3 q = p(q)d 3 q, 

i.e., 6 = J -1 — 1. By taking the divergence of the equation 
of motion (12), one obtains the equivalent of the Eulerian 
equation (3): 

,J(T,q)V x x = p(T)[J(T,q)-l]. (13) 

Of course the addition of any divergence-free displacement 
field to a solution of the previous equation will also be a 
solution. In the following, we remove this indeterminacy 
by restricting our attention to potential movements, which 
must satisfy 

V x x x = 0. (14) 

The main reason to restrict to that case is that vortical 
perturbations linearly decay, a consequence of the conser- 
vation of angular momentum in an expanding universe. 
Thus one might consider that the solutions will apply any- 
way, even if vorticity is initially present, because at later 
times it will decay away. In the same spirit, we shall mainly 
focus on growing mode solution (see Buchert and Ehlers 
(1993) and Buchert (1993) for the cases of rotational per- 
turbations and the effect of decaying modes). 

The final equation to solve obtains by rewriting the di- 
vergence of the acceleration T = x explicitly as a function 
of q 

V x T = ,J(T,q)- 1 J2^ijA ji , (15) 

where the Aij are the cofactors of the Jacobian, and the 
partial derivatives denoted by latin letter are taken with 
respect to the q coordinate (e.g., Tij = f^)- 

As in the Eulerian case, perturbative solutions are ob- 
tained by means of an iterative procedure. But this time, 
the expansion concerns the particles displacement field it- 
self, and we write it as 

# = s# (1) + s 2 * (2) + s 3 & (3) + 0(e 4 ). (16) 

The determinant of the jacobian is then similarly ex- 
panded as 

J = 1 +ejW + e 2 J^ + e 3 J^ + 0(e 4 ) 

= I +£K W + £ 2( K (2) + LiV) (it) 

+£ 3 (A K3) +jL (3) + M (3) ) + 0(£4) 
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where K^ m \ L^ m \ and M^" 1 ^ denote the m — th order part 
of the (invariant) scalars 



M = V = det[&i tj ] . 



(18) 



In other words, 

■ K (m) _ y . ^(m) _ 
2 Ei#J 



L( 2 ) = iV. J .(* a) *^-$ a .M 1) ) 

^ 2 i,i 3,3 1,3 3,1 J 

i-^l^3 V l,l 3,3 1,3 3,1 ) 

,M( 3 ) = det[<^ 



(19) 



(l)l 



The equation to solve then follows by replacing the expan- 
sions (16-18) of the displacement field and of the Jacobian 
in equations (13) and (15). We now proceed to an iterative 
solution of this equation. 

2.1. First Order 

First, keeping only the terms of order e, one gets 

A'W - = 0. (20) 

Thus K {1 \T,q) can be factorised into a temporal and spa- 
tial part, 

K( 1 \T, q ) = g 1 (T)K( 1 \n, q ), 

where K^(ri,q) is the divergence of the initial displace- 
ment field, <7i(t;) is assumed to be unity, and 

gi = k a gia + hgn- (21) 

For SI = 1, the linear growth rate g\ is given by 

gia = T~ 2 , 9lb = T 3 . (22) 

When SI < 1, 

g la = 1 + 3(r 2 - 1)(1 + rS), g lb = t(t 2 - 1), (23) 

1 T - 1 

S = - In . (24) 



2 t + 1 



The closed case of course obtains by the transformation 
(r), A, B, R, x, t, f3) — ► (irj, —A, iB, iR, —x, it, — /?). 

For a potential movement, one thus recovers Zeldovich 
solution (1970) 

* (1) (r,g) = ffi(r)«f (1) (g), 
if we define (for all m) 

Also, note that 
^\t) (X 9i {t) 



(26) 



(with initially sj 1 ^ = — VS' (</)); i- e -> the Eulerian linear 
behavior is recovered, since g\ = D [cf. Eqs. (4), (6)]. 
The logarithmic derivative of the growth factor, 



f l = a dgl 
~~ gi da 



(27) 



is useful to describe comoving peculiar velocities and 
therefore redshift distortion (see below). Near SI = 1, a 
limited expansion of the solution (24) shows that /i x 
SI 4 / 7 where x means "behaves asymptotically as". A bet- 
ter analytical fit for /i in the range 0.1 < SI < 1 is given by 
/i pa SI 3 / 5 , as was originally proposed by Peebles (1976). 

2.2. Second Order 

Collecting terms of 0(s 2 ) in equations (13) and (15), we 
find 

K^-p K ^ = -PglL^(n, q ). (28) 
Thus, the solution is again separable 

A'( 2 )(r,g) = g 2 (T)K( 2 \n, q ) (29) 
where <72(tj) = 1- The spatial part is given by 
KW(T i ,q) = LW(T i ,q), 
and the growth factor is 

g-2 = k 2 a g 2a + 2k a k h g 3h + k 2 g 2c + Lgia + hgn ■ (30) 
For SI = 1, we find 

3 3 1 

g-la = --j9l a , 92b = -^9la9lb, 9lc = ~^3lb- (31) 

When SI < 1, 

gia = l~\{r 2 -l) [r+(r 2 -l)5] 2 , 



921, = -~ A (r 2 -l) [r 3 + (T 2 -l)S] 



(32) 



g2c 



The closed case, once again, obtains by the transformation 
(rj, A, B, R, x, t, f3) — ► (irj, —A, iB, iR, —x, it, — /?). 

If only a growing mode is initially present (fcj = = 
lb), which we assume, one must have 921 9\ —> when r — ► 
00, i.e., when SI — ► 1. A limited expansion of equation (32) 
(25) then shows that l a must be equal to — |, and this yields 
the physically relevant solution 

g-2 = ~\ - \{r 2 - 1) { 1 + rS + \ [t + (t 2 - 1)5] 2 } , (33) 



which behaves near SI = 1 as 



n ~ 0" 2 / 63 „ 2 

g2 ^ -- " 9i, 



(34) 
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Fig. 2. Growth factors for Friedman-Lemaitre models with arbitrary Q and A = 0. The upper panels display the first order 
growth rate divided by the expansion factor (left), and the second order one (right), while the lower panels show the two third 
order growth factors. The solid curves show the exact solutions and the dashed curves correspond to the asymptotic behaviors 
near = 1. 



while g'j/gf — 1/2 when £1 approaches 0. 

If SI > 1, the solutions are obtained by the usual trans- 
formation, and in particular, the solution corresponding 
to (33) is 

92 = -i + ^(r 2 + l){l-r5-i[r-(r 2 + l)5] 2 }, 
S = Arctg-, (35) 



with cj2 approaching 4 — (37r/4) 2 (while g\ x 2) when SI 
goes to infinity (i.e., r — *■ 0), and with the same behav- 
ior (34) near SI = 1. 

Figure 2 shows — 3^2/ i?i versus SI, and it confirms that 
indeed g'j/gf is very nearly constant, as could be foreseen 
from the asymptotical behaviors. It is thus not surprising 
that the remaining small variation of g^ may be described 
quite accurately by Eq. (34) for 0.1 ^ SI ^ 2, which should 
cover most astrophysical uses (Bouchet et al. 1992). 




Fig. 3. Logarithmic derivatives of the growth factors of Fig. 2, i.e. for Friedman-Lemaitre models with arbitrary Q and A = 0. 
The upper panels correspond to the first order (left) and second order (right). The lower panels correspond to the two third 
order factors. The solid curves show the exact solutions, the dotted lines display their slope in Q = 1, and the dashed ones a 
better fit for 0.1 <> Q <> 1. 



The logarithmic derivative of the second order growth 2.3. Third Order 



rate, 



72 = ~r- 

92 da 



behaves as/2 x 2 SI 5 / 9 near SI = 1. Alternatively, it can be 
somewhat better approximated by 2 SI 4 / 7 for 0.1 ^£1^1 
(see Fig. 3). 



The 0(s 3 ) terms obey the equation 



^(3)_^(3) = -2^? M (3) (r . )+ ^3 (1 _i« )i (3) (7 .. )j( 3 6) 

which is separable only when SI = 1. In that case, the 
fastest growing solution is simply 



tf( 3 )(r,g) = -| fl ?(r) 



M (3 \ri, q) — j L (3 \ri, q) 
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Otherwise the third order term may be written as: 

K (3 \r, q) = g 3a (T)M( 3 \n, q ) + g 3b (T)L( 3 \n, q ), (37) 

assuming as before that g 3a {Ti) = g 3 b{ T i) = 1- Figure 2 
shows these growth rates, g 3a and g 3b , which we obtained 
by numerical integration of equation (36). Alternatively, 
by expanding the equations governing the ratios g 3a / 'gf 
and g 3b /g 3 , we find that the fastest growing parts behave 
near SI = 1 as 



21 



(38) 



Fig. 2 shows that these asymptotic behaviors near SI = 
1 provide very good fits to g 3 /gf, as was also the case 
before for the ratio g 2 jg\ [Eq. (34)]. In both cases the SI 
dependence of the ratios is very weak. 

The logarithmic derivatives of the growth factors /; = 
(a/gi)dgi/da behave near SI = 1 as 

/ 3a x3Sl 128 / 231 , /36X3S1 58 / 1 



05 



or can be better approximated for 0.1 ^ SI ^ 1 by 3 SI 23 / 40 
and 3 SI 4 / 7 respectively (see Fig. 3). 

2.4. Non-zero Lambda 

We now turn to an inflation motivated case, i.e. a flat 
Universe with a non-zero cosmological constant, 



SI- 



A 



1, 



(39) 



The expansion factor and the Hubble parameter write 



a(t) = ai sinh 2/3 ( '^H^, t 



H(t) = ifoocoth ( -H^t 



(40) 
(41) 



where a\ is the expansion factor such that SI = 1/2 and 
if 00 = ^/A/3 = H(t = 00). It is convenient to introduce 
the inflexion point For a smaller than a e 

the model-universe is dominated by matter with a neg- 
ative acceleration, whereas for larger a, it's evolution is 
dominated by the A term, the acceleration is positive 
and continuously increases. Introducing the new variable 
h = H(i)/ Hoo, the equation of motion for a particle with 
comoving coordinates x now reads 



Vj3(/i 2 -l)— + 2h — )x 



dh 2 



dh 



-28, 



(42) 



with dx/dh = dx/dh\ q . Applying the same techniques as 
in the previous part, we find the equations satisfied by the 
growth rates up to the third order 



'3(/i 2 -l)#i +2h gi 

3(h 2 -l)g 2 +2hg 2 

3(/i 2 - l)ha +2hg 3a 

Hh 2 -l)g 3b +2hg 3b 



2gi 
2g 2 - 
2g 3a 

2g 3b + 2g 3 1+ (l 



c 2g 3 l+ 



(43) 



where dots denote derivatives with respect to h. Of course, 
the spatial parts remain given by Eqs. (29), (37). 

Solutions of the first order equation are linear combi- 
nations of the respectively growing and decreasing modes 



gi+(h) = h 



dh' 



h h' 2 (h' 2 - I) 1 / 3 ' 



and gi-{h) = h. (44) 



By expressing A as a function of x = a/a e , one can verify 
that these solutions are respectively identical with D\ and 
D 2 given in LSS (eq. 13.6). The growing mode solution is 
an ugly mixture of hypergeometric functions, but we can 
identify its asymptotic behaviors: 



<7i+(a;) oc x, i<1; 



gi+{x) 



1 M 2( M 2_ 1 )l/3' 



(45) 
(46) 



We have solved the differential equations (43) recursively 
by numerical integration, to find the fastest growing solu- 
tions whose behaviors are shown in figure 4. 

Note that near SI = 1 (i.e., h — ► 00) the second order 
solution behaves as 



g2 



(SI) x -3/7 SI -1 / 143 # 2 



Similarly, the two third order growth factors behave near 
SI = 1 as 



-lv- 4/275 gl 
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J8.x-- Si g h , 3i x^ST 269 / 17875 , 3 . (47) 

Figure 4 shows that these asymptotic solutions near SI = 1 
yield good approximations for a large range of values of 
SI, 0.1 <> SI <> 1. 

The logarithmic derivatives of the growth factors (/; = 
(a/gi)dgi/da) behave near SI = 1 as 



/1 x Sl 6/n , f 2 x 2S1 153 / 286 , 
/ 3a x3Sl 146 / 275 , f 3b x 3 SI 9481 / 17875 . 



(48) 
(49) 



However, a better analytical fit for 0.1 ^ SI <> 1 would be 
(see Fig. 5) 



/1 ~ ^ 5/9 , h « 2 SI 6 / 11 , 
/ 3a «3Sl 13 / 24 , f 3b « 3S1 13 / 24 



(50) 
(51) 



Values of /1 for non vanishing cosmological constant 
have already been calculated for a flat Universe (Peebles 
1984) or in the general case (Lahav et al. 1991, Martel 
1991). Lahav et al. (1991) proposed the fit 



/1 «Sl 06 + A/70(l + Sl/2) 



(52) 



with A = A/(3ff 2 ), whereas Martel (1991) proposed the 
simpler fit 



/1 w SI 6 + A/30. 



(53) 



The computations and fits proposed by these authors are 
compared to our results in Fig. 5. 





Fig. 4. Same as Fig. 2, but for a fiat Friedman-Lemaitre models with A / (the upper left panel shows the first order growth 
divided by the expansion factor, the upper right one shows the second order growth rate, and the lower panels are devoted 
to the two third order growth factors. The solid curves show the exact solutions, and the dashed lines display the asymptotic 
behavior in Q = 1. 



3. Link to statistics if the volume V is large enough that it can be considered 

a fair sample. Then 

In order to make the connection with measurable quanti- 
ties, we proceed as follows. Let Q(x) be an Eulerian phys- 
ical quantity of interest. As usual, we shall assume that 

ensemble averages, denoted by (Q) are equivalent to aver- (Q(x)) x = (Q(q) J(q)) q , (54) 
ages over space, denoted by Q, 



Q(x) = ^ f d 3 xQ(x), 
V Jv 



where the subscripts x or q tell us whether the average is 
to be taken in Eulerian or Lagrangian space. 
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Fig. 5. Logarithmic derivatives of the growth factors for a fiat Universe with A / 0. The upper panels correspond to the first 
order (left) and second order (right) rates. The lower panels correspond to the two third order factors. The solid curves show 
the exact solutions, the dotted lines display their slope in Q = 1 [see Eq. (49)], and the dashed ones a better fit for 0.1 < Q < 1 
[see Eq. (51)]. We also show for comparison in the upper left plot the values for fi provided by Peebles (triangles), Lahav et 
al. (long dashes), and Martel (dashes-dots). 



3.1. Skewness factor in real space 

Now that large galaxy samples are becoming available, it 
becomes of particular interest to predict the moments of 
the density contrast distribution 



<^(x)) x = ((j- 1 - ir) x = ((j- 1 - ir j> q . 



The variance and skewness are then given up to the fourth 
order by 

(S 2 ) = s 2 ( JW 2 ) + s 3 (2JW J( 2 ) - JW 3 ) + 

£ 4 ( J(l) 4 _ 3 J( l) 2 j(2) + J(2) 2 + 2 J(1) J(3)) + 0( £ 5 ) 

(^ 3 ) = -e 3 (JW 3 ) + e 4 (2J(!) 4 - 3J(!) 2 j( 2 )) + 0(e 5 ), 

where all averages on the displacement field are taken with 
respect to the Lagrangian unperturbed coordinate q. 
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We consider the case of an initially non Gaussian den- 
sity field Si = eh^, but let us require for simplicity that 

the three components of the displacement field ^ are 
independent, with the same statistical law to insure ho- 
mogeneity and isotropy. We note a 2 the variance of any 
component i of the (linear) gradient field 



and S its third moment S 
duced forth moment K = < \?. 



\?;V 3 ), and K its re- 



(1)4 



3cr 4 . We thus ha 



(J(i) 3 ) = 35, ( JW 4 ) = 3K + 27a 4 = 3K + 3 ( jW 2 ) Z . 
The other term in (i5 3 ) involves the product 2 j( 2 ) 
which can readily be estimated since [after (18), (29) and 
(19)] we have 

It follows by development that (j( 1 ) 2 j( 2 )) = 6(1 + 
fi , 2/fi , i)c 4 - We thus have 



5 3 



-3S + egMK + 3(2-g 2 /gl)a 4 ) 
s gi a 4 (3-3sS/a 2 ) 2 



(55) 



S + 4-2 g 2 /gl + 2 K °\ f + 0( £ ). (56) 



3<T 4 fifi£ 



3<7 6 



For an initially gaussian field with S 1 = K = 0, we get the 
simple result 



S 3 



2g 2 /g 2 + 0(s 2 ), 



(57) 



whose first term corresponds to the pure Zel'dovich ap- 
proximation and had been found by Grinstein and Wise 
(1987). The value of the ratio of growing modes may be 
obtained from the exact solutions given before (or read 
off Fig. 2), but it is hard to imagine practical cases when 
our approximation (34) might not be sufficient. We thus 
have the handy and quite accurate formula for a gaussian 
initial field with 0.1 <> ft <> 2 (Bouchet et al. 1992) 



S 3 



28 + 6 ft~ 2 / 63 
7 ' 



which generalizes the S3 = 34/7 found by Peebles (1976) 
in the ft = 1 case. The more general formula (56) gives 
the time-evolution of the skewness factor S3 in spatially 
flat models for the class of non-gaussian initial condi- 
tions which may be generated by independent displace- 
ment fields along three axes. The general formula for the 
ft = 1 case can be found in Fry and Scherrer (1993). 

3.2. Skewness factor in redshift space 

We now turn to a generalization of the previous result (56) 
in redshift space. In redshift space, the appearance of 



structures is distorted by peculiar velocities. At "small" 
scales, this leads to the "finger of god" effect: the clusters 
are elongated along the line-of-sight due to their internal 
velocity dispersion. This is an intrinsically non-linear ef- 
fect, and we shall not be concerned with it. At "large" 
scales, the effect is reversed: the coherent inflow leads to 
a density contrast increase parallel to the line-of sight. In- 
deed, foreground galaxies appear further than they are, 
while those in the back look closer, both being apparently 
closer to the accreting structure (Sargent & Turner 1977; 
LSS, §76; Kaiser 1987). 

Kaiser (1987) estimated this redshift effect, in the large 
sample limit, on the direction averaged correlation func- 
tion, and found (6 2 ) z = (1 + 2/3 ft 6 + 1/5 ft 1 - 2 ) (6 2 ), 
where the superscript z correspond to a redshift space 
measurements. For ft = 1, {6 2 ) Z = 28/15 (<5 2 ) 

Let us thus consider the case of spherical coordinates, 
when distances to the observer would be estimated by 
means of redshift measurements. And let us now denote 
redshift space measurements by the superscript z. The 
redshift space comoving position x z of a particle located 
in r(q) = ax(q) is x z = r/(aH) (with H = a/a). The 
real space perturbative expansion (7) is then replaced by 



+ [l + /i(*)]</i(*)£ (1) (g) 

+ [l + h(t)]92(t)^ 2 \q) + 0(e s ) 



(58) 



where we have explicitely used the separability of 3''- 1 -' = 

gi{t)& (1 \q) and «F (2) = g 2 (t)^ (2 \q) [eqs.(2.1) and (29)]. 
In the limit of an infinitely remote observer, say along 
the r3-axis, the observed density constrast 6 Z in comoving 
coordinates is simply 6 Z (xi, x 2 , x z 3 ), which amounts to ap- 
proximate spherical coordinates by cartesian ones. All we 
have to do, then, is to replace everywhere in the calcula- 
tion of S3 



by 



(m) 



9m ^ 



(l + f m )g m &3 m) 



for m = 1 and 2. We have, 

(jW 4 )= [2 + (l + / 1 ) 4 ] A' + 3(j( 1 ) 2 ) 5 



with 



(jW 2 )= [2 + (l + / 1 ) 2 ] a 2 



This shows that in our so-called "infinite observer limit" , 
we have (8 2 ) z = (1 + 2/3 /1 + 1/3 f 2 ) (8 2 ) which slightly 
differs from Kaiser's calculation (who did a calculation 
in spherical coordinates instead of our rectangular ones). 
Indeed, we find that the redshift space variance is boosted 
by a factor 30/15 for ft = 1, instead of his value of 28/15. 



14 



F.R. Bouchet et al.: Perturbative Lagrangian approach to gravitational instability 



The term (jW 2 ^ 2 )) is now equal to 
[2 + 4(1 + /i) 2 + (g 2 /gl) (2 + 4(1 + ft) + 2/ 2 + hh U)] 



with 



^(1)^(1)^(2) 
*1, 1*3, 3*3, 3 



(2 + /i)(*. 



-,(1)2^(2) 



containing all symmetry breaking terms. To evaluate U, 
we use discrete Fourier transforms, which we denote by 
hats, in a very large volume. It follows from (28) that 



= (w^!/ifci 2 )EE 

k 1 i>j 

[ki(k J -k')-k'(k l -m^ i \k , )^ i \k- 



k') 



SI = 1, we have = 1 = / 2 /2 (and g 2 lg\ = — 3/7), which 
yields for Gausssian initial conditions S3 = (35 + £)/7 
,while, for SI = 0.1, S 3 « (34.5 + 0A£)/7 (to be compared 
with the value of 34/7 in real space). 

The formula (59) obtained above applies only in the 
limit of large volumes (like Kaiser's result), since we have 
taken the limit of an infinitely remote observer. A full 
calculation along the previous lines, but in spherical co- 
ordinates (and including the effect of smoothing) may be 
found in Hivon et al. (1994). In any case, equation (59) 
clearly shows that the ratio S3 is, for gaussian initial con- 
ditions, nearly independent of the value of £1, nor is it 
affected by redshift space distortions. It is interesting to 
note, though, that in the non-gaussian case, the distortion 
might be rather large, for large enough S or K. 



(with = *$(fc)exp(ifc • x)). S 



(2)/ 



^(fc)*^')) = 8 K (i - j)8(k - k')P(k), 



i.e., is zero for i ^ j (8k is the Kronecker symbol), or 
k ^ k', it follows that (^^V*^*^) is zero for i = 3 
and is equal to 

(92/9l) E^3 2 (^3 - k' 3 f/(k - k'fP(k)P(k') 

k,k' 

otherwise. Since g 2 < (at late times when the growing 
mode is dominant), we have 



similarly, 



4. Comparison with other approximations for 
spherically symmetric perturbations 

So far, we have mainly considered rigorous uses of the La- 
grangian perturbative approach, for instance the deriva- 
tion of a second-order quantity, the skewness of a PDF 
with the help of second order perturbation theory. Now, 
we examine to what extent the second-order theory brings 
improvement to Zeldovich approximation, when both are 
used as approximations to the real dynamics, i.e. outside 
of their rigorous validity range. In this section, we first ex- 
plicitly derive Eulerian expressions from their Lagrangian 
counterparts (up to second order). Then we use the spher- 
ically symmetrical model to compare various perturbative 
approaches with the exact solution. In the next section, 
we shall compare to numerical simulations from Gaussian 
initial condition with a power-law power spectrum. 



^(1W1W2)\ < Q 

*i, 1*3, 3*2, 2/ ^ u ' 



which implies 

sMMK!) > MM3W 2 )) = (g,/gi)a\ 



As a result, U is bounded, negative, and lies between 
and g'j/gl- It follows that 



6 



(2 + (l + / 1 ) 3 ) S 
£ (2 + (l + / 1 )2) 2 ,74 

2 (2 + (l + / 1 ) 3 ) 2 5 2 

(2 + (l + / 1 )2) 3 
2 (2 + (l + / 1 ) 4 ) K 

(2 + (l + / 1 )2) 2 ,74 
e 1 + 2(1 + /1) 2 + ( g2 /ffi) P + 2/1 + h + AML x 
[2 + (l + /i) 2 ] 2 



with 1 > £ > 0. Of course, we recover the real space 
result (57) if we set fi = f 2 = 0. On the other hand, if 



4-1. Lagrangian/ Eulerian approach up to second order 

Knowing the fastest growing solution up to second order 
in the Lagrangian formalism, we are now able to find its 
Eulerian counterpart. In this section, it is more convenient 
to slightly change the previous notations: we drop the e (or 
equivalently the displacement field is redefined to include 
the epsilons) and we use gj i = gj(r = r 8 ) = gj(t = t{), 
j = 1,2 (i.e. we do not impose gj(t = t{) = 1 anymore). 
The Eulerian initial conditions are generally specified by 
an initial density profile 



8 i (x) = 8(x,t i ) = 0(e), 



(60) 



and an initial velocity profile. We assume, as before, that 
Si is dominated by the fastest growing mode. Let us also 
recall that, at lowest order in e, the initial displacement is 
determined by 



V, • ¥(ti,q) = V, • * (1> (q) = -6i(xi[q]), 
V 9 X# (1) (g)=0. 



(61) 
(62) 
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The mass conservation equation (8 = J 1 — 1) yields 



6(x[q}) 



9i 

9i,i 



9l 



ll 
9h 



■ ff f(v,* (1) ) 2 

1,8 

.¥.(2) 



2l 

9l 



where (q) is the curl-free quantity verifying 



(63) 



mode of the Lagrangian perturbative expansion (7). What 
we will call Eulerian linear theory corresponds to the first 
term in right hand of Eqs. (68), (71). 

Strictly speaking, the second order expressions given 
above are valid only if 



9\l9\,i > 1, 
£9i/gi,i < 1- 



(72) 
(73) 



2 ^ 



IJ m.m 



l m m.1 



(64) 



It is then useful to introduce the Eulerian displacement 

~ ( le ) 

(a;), defined by 



* ile \ X ) 



x). 



(65) 



By construction, this displacement is curl-free in Eulerian 
0. It is simply related to the initial 
-8i(x). We can write 



~(le) 



space, i.e. V x 

~ (le) 

density contrast by 



v,.* (1) 



# (1 %)-^V„# (le \# (le) + 0( £ 2 ), (66) 

9l,i 



v,.* (le) - 



9i 

9i, i 



v x (v x .* (le) 



>} 



The Eulerian density contrast 6(x) then is: 
«(x) 



The first condition insures that transient modes and sub- 
dominant modes are negligible. The second condition is 
necessary for the perturbative solution to be valid. How- 
ever, as we will see in the next paragraph when we consider 
the spherically symmetrical model, the condition (73) may 
be relaxed in the Lagrangian case, as far as the density 
contrast is concerned. Indeed, mass conservation is by con- 
struction verified in the Lagrangian approach. It explains 
why Zeldovich approximation, although a first order (lin- 
ear) solution, provides a very good qualitative description 
of the density distribution in pancake models. In the Eule- 
rian case, though, the condition (73) has to be obeyed, oth- 
erwise the approximation (68) becomes completely wrong. 
With regard to velocity fields, one can expect the Eule- 
rian perturbative approach to be almost as accurate as 
< ~ le ' ) -|-0(£ 2 ).(67)he Lagrangian one at the same order, since mass conser- 
vation is less important. For example, roughly speaking, 
expressions (70) and (71) are equivalent when taken at 
first order, if (70) is evaluated at x' = q. 



9i, 

il 

ah 

ll 

ah 



[(v*.# 



,f.( le ))2 



{v,(v,.* (le) )}.* (le) ] 
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9\ 



v,.* (2e) + o( £ 3 ). 



£( 2e ) 



The quantity 

1 \ - 



(a;) is curl-free and verifies 



v,.* (2e) 



2^ 



~(le)~(le) 
1,1 m,m 



~(le)~(le) 

-■- / t-t-i -■- t-t-i l 



(68) 



(69) 



A similar transformation can be made for the velocity 
field. In the Lagrangian approach, the peculiar velocity of 
an element of matter is written as 



v(g,*) 



ll 

9i, 



-* (1 V 



11 

at 



0(e 3 ). 



(70) 



This leads in Eulerian coordinates to 
v(x,t) 



11 

91,; 

11 

ah 



■* (le \x) 



■* (2e \x) 



2im x ^\^\ 



0(e 3 pi) 



For further reference, let us recall that Zeldovich ap- 
proximation is nothing else than the fastest linear growing 



4-2. Example: the spherical model 

The spherical model can be exactly solved analytically 
before shell-crossing. It thus provides a good test for var- 
ious perturbative models. The aim of this paragraph is 
to compare different approximations with the exact solu- 
tion. Firstly, we consider the case of the top hat model, 
which is equivalent to look at the center of the perturba- 
tion. The density contrast will be studied, as well as the 
divergence of the velocity field. Then, the density contrast 
and the velocity field of a given profile will be analyzed 
as functions of the radial coordinate r. The perturbative 
solutions studied here are the standard Eulerian linear ap- 
proximation, the Zeldovich approximation, the Eulerian 
Second order approximation and the Lagrangian second 
order approximation. Details of the solutions for the evo- 
lution in the various approximations are given in appendix 
A. In Section 4.2.1, we also consider Bernardeau's model 
(1992). This model provides an exact relationship between 
the density contrast and the divergence of the velocity field 
in the weakly nonlinear limit (when the variance (6) 2 of 
the distribution is infinitely small). It is based on a sta- 
tistical approach, so it may be somewhat meaningless to 
use it here, since we consider a particular symmetry. How- 
ever, it will be interesting to notice that the calculation 
of Bernardeau provides a very good fit of the relation be- 
tween V.v and 6 for the spherical top hat model, which 
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Fig. 6. Computed density contrast 8 as a function of the exact solution 8e for various approximations for a spherical top hat 
perturbation (Q = 0.1). The left panel and the right panel respectively correspond to 8e < and 8e > 0. The abbreviated 
labels in the right panel correspond to the following models: "Ex." stands for the exact solution, "El" for Eulerian linear 
theory, "E2" for Eulerian second order perturbation theory, "LI" for Zeldovich approximation, "L2" for Lagrangian second 
order perturbation theory, and "B" for Bernardeau's model. 



is rather useful when one wants to use it as a toy model, 
since exact solutions have rather complicated expressions. 

In the following, we assume that shell-crossing has not 
taken place. Results are given for SI = 0.1 (and A = 0), but 
the qualitative behavior of the solutions does not change 
significantly with £1 and the conclusions given here can be 
generalized for any reasonable value of SI (0.05 ^ SI ^ 3). 

4.2.1. The top hat model 

Let us consider an initial density profile given by 



8(x^,t^ — Si , Xi ^ x\^, 
f>(%i,ti) = 0, Xi > xi ti . 



(74) 2 



According to Birkoff's theorem, the evolution of a La- 
grangian sphere of initial radius Xi < x\ y i depends only 
on its contents. Then the matter distribution inside the 
sphere of radius xi(i), with x\(ti) = x\ y i remains homo- 
geneous while the system is evolving. For x < xi(i), the 
density contrast 8(x,i) is thus a function of time and 6{. 
This reasoning can be generalized to any initial density 
profile 8(xi,ti), when one considers the center of the per- 
turbation: 8(0, t) depends only on 8(0, ti) and t. 

Figure 6 gives the density contrast 8 as a function of 
the exact solution 8e for various perturbative models and 
for Bernardeau's model. As expected, the Lagrangian ap- 
proach seems to be much more efficient than the Eule- 
rian one. Even Zeldovich approximation, which is only 
valid at first order, is better than the Eulerian second 
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Fig. 7. The quantity — V x .v/(aH8) as a function of 8, in loga- 
rithmic coordinates, for various approximations in the spherical 
top hat model (Q = 0.1). Notations are the same as those used 
in Fig. 6. The Lagrangian second order approach (L2) and the 
Eulerian second order approach (E2) cannot produce arbitrar- 
ily small density [see Fig. 6, Eqs. (75), (76)], which explains 
the discontinued line for E2 and the vertical asymptote for L2. 




Fig. 8. The density contrast (left panels) and the velocity field (right panels) versus the distance to the center of the halo for 
two initial values of 8 (see text). The top and bottom panels respectively correspond to 8 = 0.002 and 8 = 0.005. For a better 
legibility, the bottom left graph is in logarithmic coordinates. 



order approximation (hereafter L2). Note that at second 
order, arbitrarily underdense regions cannot be obtained. 
In the Eulerian case, the largest underdensity that can be 
reached is 



,E2 



28(14 + 3fi- 2 / 63 ) 



-0.3. 



(75) 



For Lagrangian second order perturbation theory, we have 
instead 



,L2 



7_ 

12' 



fi 2/63 



-0.75, 



(76) 



which is of course much better. In the regime —0.8 iS 6 iS 
3, the accuracy of the Lagrangian second order approxima- 
tion is better than 7%. Although it considers a statistical 
average and not a single object, Bernardeau's model is 
quite good. 





Fig. 9. Same as in Fig. 8, but for a void. Even for a moderate final density contrast (8 ~ —0.4), the Eulerian second order 
approach gives bad results. 



Figure 7 is similar to Fig. 6, but it gives the quantity 
Vx.v/(aHi) as a function of 6. As in Fig. 6, the model 
of Bernardeau gives a very good approximation of the ex- 
act solution. One can see that a second order solution 
is needed to get the appropriate slope of V x .v /(aHb) in 
the vicinity of 6 = 0, as expected. Indeed, this slope corre- 
sponds to the first correction to the horizontal dotted line, 
which refers to the Eulerian linear theory. In this case, Zel- 
dovich approximation is not really better than the Eule- 
rian linear theory. This is not surprising, as already argued 
in Sect. 4.1. Similar conclusions can be reached when one 



considers second order approximations. The discontinued 
line in the Eulerian second order case and the vertical 
asymptote in the Lagrangian second order case both re- 
flect the existence of a maximum underdensity that can 
be reached in these approximations. 

4.2.2. Study of a given profile 

The spherical top hat model is quite extreme. On can in- 
deed expect that perturbative models fail in approaching 
the exact solution at the center of the fluctuation. As we 
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shall see through an example, for a given initial profile 
6i(ri), the Lagrangian second order approximation is able 
to describe density contrasts as large as ten, if one consid- 
ers regions far enough from the center of the fluctuation. 

Following Martel & Freudling (1991, hereafter MF), 
we take as initial profile 



1-th 



n, 



£.n, 



(77) 



The typical radius of the initial perturbation is thus given 
by n,i, while £ parameterizes the width of the transition 
between the regime 6;(r;) ~ 8 and the regime 6;(r;) ~ 0. 
We choose, as MF, 



n, 



1, £ = 0.3, a/cii = 1500. 



(78) 



Figures 8 and 9 give for a spherical halo and a void, 
the density contrast and the velocity field as functions 
of r/ri(t). The function r\(t) is the Lagrangian position 
of points initially at r\i, or in other words, the typical 
radius of the fluctuation at t. The peculiar velocity v is 
given in units of the Hubble flow for r = ri(t). We take 



|<5i| = 0.002, \6 2 \ = 0.005, 



(79) 



to test the various approximations: \6i\ is such that the 
final value of 6 is moderate, about unity for a halo and 
about —0.4 for a void; | <5 2 1 leads to large density contrasts, 
about 40 for a halo and about —0.6 for a void. The calcula- 
tions are assuming Slo = 0.1, but the results are quite sim- 
ilar for different Slo, although one has to decrease the value 
of 6 when Slo increases, in order to get the same final 6. The 
values of Slo and a/a; choosen here give gi/gi t i — 290 1, 
so the validity condition (72) is verified. But condition (73) 
is not verified, since we get £g\/g\ t i — 0.58 for \8\ = 
and sgi/gi t i ~ 1.45 for \6\ = l^l- This is outside the ap- 
propriate regime for Eulerian second order perturbation 
theory at least for \6\ = | <5 2 1 - 

Figs. 8 and 9 confirm the results of the previous para- 
graph, in terms of relative accuracy of the various approx- 
imations for moderate final density contrasts: 

density contrast : 

Lagrangian second order > Zeldovich 2; (80) 
Eulerian second order > Eulerian linear theory, 



Velocity field : 

Lagrangian second order 
Zeldovich 



Eulerian second order > (81) 
Eulerian linear theory. 



For large final density contrasts, the Eulerian approach 
becomes particularly inefficient, while the advantage of 
the Lagrangian approach is much less pronounced for the 
velocity field. 

Overall, the second order Lagrangian approach gives, 
for moderate final 6, an excellent approximation of the 
density contrast and the velocity field. It appears to 



correctly reproduce density contrasts as large as ten. 
Thus these comparisons suggest that by releasing con- 
dition (73), i.e. by extrapolating the perturbative La- 
grangian solutions and using them as approximations to 
the real dynamics, one gets a good idea of the behavior of 
the system in the non-linear regime. 



5. Comparisons with a Numerical simulation with 
Gaussian initial perturbations with a power 
spectrum oc k~ 2 

In order to see whether the conclusions reached in the 
spherical case hold in a more realistic setup, we now study 
the evolution of three-dimensional, Gaussian, initial per- 
turbations with a power spectrum oc k~ 2 . A similar study 
(with similar conclusions), but for truncated power spec- 
tra may be found in Melott et al. (1994). We compare the 
evolution of 64 3 particles in a 128 3 PM numerical sim- 
ulation with the evolution computed by using the first 
(Zeldovich, or LI) and second order (L2) Lagrangian ap- 
proximations. Periodic boundary conditions are assumed. 
All Eulerian quantities such as the density and the velocity 
fields are computed on a 64 3 grid. 



5.1. A Slice through the data 

The figure 10 shows the particles positions in the same 
thin slice of the box, Li, ox /64: thick, when the simulation 
(left column), the LI (Zeldovich, middle column), and L2 
(right column) approximations are used to compute the 
evolution of the same initial conditions (hereafter IC). The 
rows from top to bottom show the positions after an ex- 
pansion by a factor a = 4, 8, and 16, respectively. 

Of course, the approximations fail to describe the inner 
structure of regions where shell-crossing has occurred, and 
it is interesting to rather compare the resulting smoothed 
density fields. The first three rows of figure 11 display the 
density contours corresponding to the slices of figure 10. 
As is already well-known, Zeldovich approximation, de- 
spite its simplicity, does very well at describing the general 
large-scale texture of the density field. We also note that 
at this purely visual level, L2 does not appear to bring 
very marked improvements over LI. 

In both approximations, the trajectories are entirely 
determined by the initial local velocity (LI) and also the 
acceleration (L2) field. Thus, contrary to the real case, 
if small scale variations are initially present, they will get 
amplified and transferred to much larger scales (even more 
so for L2, see figure 1), since the recall force of a forming 
non-linear clump is absent. Kofman et al. (1992) first pro- 
posed to truncate the initial power spectrum of its small 
scale power to improve the results of Zeldovich approx- 
imation. And indeed Coles et al. (1993) and Melott et 
al. (1993) showed that this smoothing of the initial con- 
ditions, at a scale close to the scale of non-linearity when 
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Fig. 10. Slices Lbox/50 thick of a model evolved from n = —2 gaussian initial conditions with respectively a PM numerical 
simulation (left column), Zeldovich approximation (middle column) and the Lagrangian second order approximation (right 
column). The rows from top to bottom correspond to a = 4, 8, and 16. 
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Fig. 11. The tree top rows show the density contours for the same slices than in the previous figure (i.e. PM, LI, and L2 
are from left to right, and a =4, 8, and 16, from top to the third row). The left figure of the bottom row show the computed 
evolution of the variance (in cubic cells of size 1/64), when the initial conditions are smoothed with Gaussian filters of various 
scale £ s . The adjacent figures show the corresponding a = 16 density contours for the LI and L2 case, when an initial smoothing 
has applied at a scale £ s = 1/64, to be compared with the unsmoothed case directly above. 
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Fig. 12. Density measured along a one-dimensional scan in a box evolved with a PM simulation (solid), and by using the LI 
and L2 approximations (dots and dashes respectively) at a=16. In both cases, the density has been obtained by smoothing the 
data with a gaussian filter with £ s = (see text). The top panel compares the results for identical initial conditions, 

while on the bottom the initial conditions for the LI and L2 approximations have been smoothed with the same £ s = Lb ox jM 
gaussian filter used at a=16. 
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Fig. 13. Approximate density contrasts, 62, obtained by LI (left) and L2 (right) versus the PM one, 6\. The top panel is for 
a = 8. . The middle and bottom panels are at a = 16 and correspond respectively to unsmoothed and smoothed IC. See text 
for details. 



24 F.R. Bouchet et al.: Perturbative Lagrangian approach to gravitational instability 




Q I i i i i I i i i i I i i i i I i i i i I 

12 3 4 

Log 10 N 

Fig. 14. The CPDF P N (£) at a = 16, for a PM, LI, and L2 "mover" (solid, dots, and dashes respectively). From left to right, 
the curves are for the scale log 10 £=-1.6, -1.4, -1.2, -1.0, -0.8. The corresponding variances (see eq.(85)) measured in the PM 
simulation are (<*> 2 ) — 3.3, 1.6, 1/1.26, 1/2.5, 1/5. Note that for a greater legibility, each curve has been offset as compared to 
its left neighbor by +0.5 on the y-axis. 




Fig. 15. The third, forth and fifth normalized moments of the density distribution for a = 8 (left) and a = 16 with smoothed 
initial conditions for approximations (right). The squares correspond to the simulation results, while triangles and circles denote 
the results obtained using respectively Zeldovich approximation and L2. 
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5 10 1 

Fig. 16. Distribution function P v (<t) at a = 16, for a PM, LI, and L2 "mover" (solid, dots, and dashes respectively). The top 
pannel shows a comparison for unsmoothed intial conditions, and the bottom one for smoothed ones (with a gaussian filter of 
size £ s = i(, ox /64).In each panel, the curves from top to bottom, are for the following <j 2 pm °f the PM simulation : 3.3, 1.6, 
1/1.26, 1/2.5, 1/5. For a greater legibility, each curve has been offset as compared to its neighbor above by -1 on the y-axis. 
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Fig. 17. Same as figure 13, but for the modulus of the velocity. LI is on the left, L2 on the right. The top row is at a=8 and is 
not smoothed. The bottom rows corresponds to a latter time, a = 16, when the measurements are performed after smoothing 
of the final velocity field, in the case of unsmoothed (middle) and smoothed (bottom) initial conditions. 
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Fig. 18. Distribution of angles between the true velocity, and 
the one computed using LI (dots), or L2 (solid). The top plot 
shows the result (for unsmoothed IC) at a = 8. The other plots 
show the result after smoothing of the final velocity field, at 
a = 16, for unsmoothed (middle) and smoothed IC (bottom). 



the comparison is made, improves dramatically the pre- 
dictions of Zeldovich approximation. 

One may think of several ways of choosing the exact 
value of the scale at which the initial data is filtered. One 
would like to choose the smallest scale that still prevents, 
at the time the comparison is made, the unphysical "blow- 
up" of small scale structures. One way to do this is to 
choose that scale such the variance of the smoothed initial 
field is maximal at that time. The left panel of the bottom 
row of figure 11 shows the time evolution of the variance 
for various initial gaussian smoothing scales. It suggests 
that at a = 8, no smoothing should be necessary. And in- 
deed we found no improvement in all the tests we consid- 
ered when we used an initial smoothing with a filter of size 
4 = L box /64 (i.e. the filter W(x) is oc exp(-x 2 /2£ s )). But 
the same figure suggests that this smoothing scale should 
be appropriate at a = 16. The adjacent middle and right 
column figures of the bottom row display for reference the 
density contours at a = 16 for the LI and L2 case, when 
this initial smoothing has been applied before computing 
the particles trajectories. In the following, we shall refer 
to these cases as LIS and L2S. 

This visual comparison suggests that indeed an initial 
smoothing leads to a better, more contrasted, description 
of structures at late stages. We now proceed to a finer 
analysis of the computed evolution. 

5.2. Density scans along a line 

The figure 12 shows the density along a line, for the dif- 
ferent approximations at a = 16. This line was chosen to 
go through the densest structure present in the simula- 
tion. The density was obtained by first affecting the par- 
ticles to a 64 3 grid using a cloud-in-cell interpolation, and 
then convolving the resulting field with a gaussian filter of 
£ s = 1/64. The top panel compares the results for identical 
initial conditions, evolved with the PM (Solid), LI (dots) 
and L2 (dashes) movers. The bottom panels compares the 
LIS and L2S density scans when the initial conditions for 
the LI and L2 approximations have been smoothed with 
the same £ s = Li, ox /64: gaussian filter used to obtain the 
final field, at a=16. 

For relatively weak density contrasts, 3 ^ p ^ 1/3, L2 
is marginally better than LI, and both are much better 
than LIS and L2S. In a strongly underdense regions (at 
Y ~ -55), LI overestimates the density contrast, while L2 
underestimates it. For smooth initial conditions, both LIS 
and L2S underestimate the real density contrast, but of 
course LIS is much better than L2S, since LI was initially 
overestimating it. In a strongly clustered region (at Y ~ 
.95), both LIS and L2S do much better than LI and L2 
(i.e. both the peak width and height are closer to the PM 
values), and L2S clearly surpasses LIS in terms of peak 
width and height. 

These scans therefore suggest that smoothing is only 
appropriate for improving results in very dense regions, for 
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which L2 appears superior to LI, the reverse being true in 
the underdense ones. These findings thus tend to confirm 
the results obtained in the spherical case (see Eq. (76)). 

5.3. Density Cross-correlations 

A cross-correlations of the density fields computed using 
LI or L2 with the simulation allow a more quantitative 
appraisal of their respective merits. Figure 13 shows the 
approximate densities, 1 + 62 versus those in the simu- 
lation, 1 + 61 at a = 8 and a = 16. In all cases, the 
density is obtained by a CIC interpolation to a 64 3 grid, 
and concentric lines surround respectively 25, 50, 75, 90 
and 99 % of the points. The top row shows the results 
at a = 8. The two other rows correspond to a = 16. 
Both have been smoothed at that time with a gaussian 
filter of £ s = Lb 0X /QA. The middle row corresponds to 
unsmoothed initial conditions, while the bottom one cor- 
responds to LIS and L2S, i.e. the initial conditions have 
been smoothed with the same gaussian filter than the one 
used at the time of comparison. 

These plots do confirm the previous findings and show 
that L2 gives better results than Zeldovich for both expan- 
sion factors, in particular for high density contrast 6 ps 10, 
and with a smaller dispersion around the true values. The 
correlation coefficients 



Ci 
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(H 



(S 2 2) 



id Co 
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provide global measures of the performance of the approx- 
imations. They are given in the upper left corner of each 
plot and comfort our conclusions. But by very definition, 
they are mostly sensitive to the densest areas, which hides 
the fact that L2 is worse than LI to describe accurately 
the lowest negative density contrasts. 

5.4. Count Probability Distribution Function 

We now turn to the statistical description of large scale 
structure in terms of the count probability distribution 
function (CPDF) and its first few moments. The CPDF 
allows a simultaneous appraisal of the approximations at 
different smoothing scales. The CPDF Pn(£) is the prob- 
ability of finding N points in a sphere of radius £. It is ob- 
tained by throwing at random a large number of spheres 
and by recording there occupation number. As we shall 
see below, the moments of the CPDF are directly related 
to the correlation functions (6®) = £q introduced in § 3. 

5.4.1. The CPDF 

Figure 14 shows the measured Pn(£) for different £ at 
a =16, both for unsmoothed (top) and smoothed (bot- 
tom) initial conditions. The PM result is shown by a solid 
line, while dots and dashes show those corresponding to LI 
and L2, respectively. The measurements scales are, from 



left to right, log 10 £/L box = -1.6, -1.4, -1.2, -1.0, -0.8. The 
corresponding variances (see eq.(85)) measured in the PM 
simulation are (<S 2 ) ~ 3.3, 1.6, 1/1.26, 1/2.5, 1/5; they 
span all the range of the transition from the weakly to the 
strongly non-linear regime. 

The figure shows that, in the unsmoothed case, L2 cor- 
rects LI in order to make the corresponding approximate 
CPDF a quasi-perfect description of the PM CPDF, apart 
from the underdense regime (which are less probable in 
L2 than in PM, i.e. we recover that L2 underestimates the 
negative density contrasts), and the densest regime where 
the improvement brought by L2 over LI still does not al- 
low to describe the highest and rarest density peaks. In 
the case of smoothed initial conditions, the Ll-CPDF and 
the L2-CPDF change according to what could be expected 
from our earlier analysis: LI now underestimates also the 
negative density contrasts; the description of weak density 
contrasts by both LI amd L2 is somewhat degraded, but 
smoothing does bring an improve description of the dens- 
est part of the simulation. It will prove convenenient to 
rephrase those statements in terms of values of the num- 
ber of standard deviation, which we now compute, as well 
as higher moments. 

5.4.2. Moments of the CPDF 

The central moments of order Q of Pn(£) are defined by 



{£) =< (N-Nf >= Y^(N -N) Q P N (£), (82) 



N=0 



with 

77 = n£ 3 . 



(83) 



This in turns allow to compute the Q-body averaged cor- 
relation functions 



Q / r! 3 r 



ri... / d rQ^ Q (ri, ...,r Q ). 



(84) 



Indeed, £ Q (£) is given for Q < 5 by (Fry & Peebles 1978, 
Peebles 1980 for the four first orders) 



Nt 2 (£) = V2-N, 

N%(£) = a* 3 - 3a<2 + 2#, 

N 4 £ 4 (£) = m- 6/13 - 3a*| + 11^2 - ON, 



(85) 
(86) 
(87) 



N ( 5 (£) = //5-10//4-10//2/"3+35// 3 +30//^-50// 2 +24#.(88) 

These expressions take the discreteness of the set into ac- 
count. The larger is Q, the larger are the contributing 
values of N in the sum (82). The high-order correlations 
are thus related to high density contrasts. 

As shown earlier in § 3, in the weakly non-linear 
regime, the term of order n — 1 in the density contrast is 
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required to compute (Cn(£)}- And therefore, the term of or- 
der n — 1 in the Lagrangian displacement field is necessary 
to make a correct estimate of the S n : LI gives exactly the 
evolution of the variance in the weakly non-linear regime, 
and a fraction of the correct answer for the skewness or 
the kurtosis of the PDF. L2 yields the correct weakly non- 
linear answer for the skewness. 

Figure 15 shows the normalized third, forth and fifth 
moment of the density distribution measured in the sim- 
ulation and in the LI and L2 cases, as functions of the 
computed variance ^2- The measurement scales I lie be- 
tween —2.2 < log (£/ Ly MX ) < —0.8. At large scales, we 
do recover what is expected form the perturbation the- 
ory. But we also see that L2 gives reasonably accurate 
values for S5 , which should require in principle a third or- 
der calculation. Even more insteresting is the transition 
regime, at scales £ such that (Cn(£)} ~ 1- Here, even if 
each approximation under-estimates the second moment 
at scales progressively more strongly non-linear, the ratios 
S3, S4, S5 remain much better estimated by L2, up to vari- 
ance 2. It was already known that LI gives, even when 
(Cn(£)} ~ 1 a good qualitative results, even outside of its 
proper validity range. Figure 15 shows that the "miracle" 
extends to L2, with the added advantage of giving a much 
better quantitative description than LI, up to S5 . In view 
of this result, it seems doubtful that third order theory L3 
would give appreciably better results than L2. 

5.4.3. Standard Distribution 

Finally, we can use the computed variance (which includes 
discreteness effects, but they are negligible for all scales 
but the smallest one shown) a = H2 in order to produce 
normalised distribution function P v (<r), with v = (N — 
N)/a. The corresponding curves are shown on figure 16. It 
shows again that smoothing helps only in the description 
of rarest events, and is otherwise rather detrimental. A 
rough measure of the quality of the approximation may be 
obtained by searching, for a given a, what is the largest v 
for which a particular accuracy criterion is met, e.g. that 
the approximate CPDF does not deviate by more than 
a factor of two from the PM one. The results of such a 
measure are given in table 1. 

5.5. Velocity Field 

We now analyze the velocity field, in terms of moduli and 
angles. Figure 17 shows the modulus of velocity in the two 
approximations versus the true velocity modulus. The ve- 
locity are in units of the Hubble velocity across the box. 
As before, at a = 16, the velocity field was smoothed be- 
fore its modulus was computed (both for unsmoothed and 
smoothed initial conditions with a gaussian filter of size 
Lbox/Q^)- In all cases, L2 provides a substantial improve- 
ment over LI, which tend to overestimate the velocities. 
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Table 1. The largest number, vm, of standard deviations, a 
for which the approximate CPDF does not deviate by more 
than a factor of two from the PM one. The first column gives 
the PM variance, <7p M ,the following ones give vm for LI and 
L2 applied to unsmoothed initial conditions, and then LIS and 
L2S, i.e. LI and L2 applied to smooth initial conditions. These 
measurements were made at a =16. second and third one give 



We can also see that smoothing of the initial conditions, 
is on this test clearly beneficial. 

Let us now define a misalignment angle descriptor, 

C (.) Sg .(V.(.),V^)) = ^ff , (89, 

where V s (x) and V a (x) are respectively the simulated 
and approximated velocity measured at the same point x. 
This function c clearly lies between -1 and 1 and we fur- 
ther define Af(c)dc as the number of points in the interval 
[c, c + dc]. Only < c < 1 is showed, but the number of 
points with negative cosine is negligible and continuously 
decreases with decreasing c. 

For two uncorrelated velocity fields Af(c) would be a 
constant independent of c, whereas positively correlated 
fields Af(c) would lead to a peack at c = 1. Figure 18 shows 
this distribution function for LI and L2. At a = 8, the 
velocity field given by L2 is a little more correlated with 
the simulation than the one given by Zeldovich, whereas 
for a = 16 the two approximations give rather similar 
results. And we find again that at a = 16, smoothing of 
the initial conditions does improve the comparison with 
the PM results. 

The improvement brought by L2 as compared to LI 
is less pronounced for this misalignment angle than for 
the velocity moduli. But overall L2 gives a much better 
estimation of the velocity field than Zeldovich approxima- 
tion. Zeldovich approximation gives good results for the 
density evolution mainly because, as any consistent La- 
grangian approach, it conserves mass and momentum at 
first order; the good L2 results shows the improvement 
brought by conserving them at second order. 

6. Conclusions 

In this paper, we extended further the perturbative La- 
grangian approach initiated by Moutarde et al. (1991). 
We gave in particular the growth rates (for the fastest 
growing modes) up to third order for the matter era of 
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Friedman-Lemaitre models with arbitrary density param- 
eter £1, with either a zero cosmological constant, A = 0, 
or in the spatially flat case, SI + A = 1. 

We used these results to compute the skewness of the 
density field, (6 s ), in the weakly non-linear regime, both 
for gaussian and a broad class of non-gaussian initial con- 
ditions. We then showed how to compute the same quan- 
tity in redshift space, when the "particle" Eulerian posi- 
tions, v are estimated by r/H. It shows that the relation 
between skewness and variance, which is unaffected by this 
transformation in the gaussian case, may be quite modi- 
fied for non-gaussian initial conditions, depending on the 
strength of the initial deviations from "gaussianity" . 

We then checked whether the second order correction 
to Zeldovich (linear) approximation improved the well 
known ability of that approximation to follow some as- 
pects of the dynamics into the non-linear regime, even 
when the validity conditions of these approximations are 
not fulfilled anymore. We found that this is indeed the 
case, both for various spherically symmetric perturbations 
of known evolution and by detailed comparisons with a 
numerical simulation starting from gaussian initial condi- 
tions of spectral index n = —2. 

The main improvements brought by the second order 
correction concern the description of regions with large 
density contrasts (8 1) and of the velocity field. We 
found that the density PDF is computed with reasonable 
accuracy up to as much as 18 standard deviations, when 
the latter reaches about 1.8. And indeed the five first or- 
der moments of the PDF are accurately estimated till that 
stage, while Zeldovich approximation would quite under- 
estimate them. The improvements brought to the descrip- 
tion of the velocity field are also rather spectacular, as may 
be judged for instance by the cross-correlation coefficients 
C\ and C'i- 

Given that the second order correction to Zeldovich 
approximation retains the advantages of the latter, i.e its 
ability to derive any stage of the evolution by simple cal- 
culations performed at an initial stage, we conclude that 
this second order Lagrangian approximation will be quite 
useful for the understanding of large scale structures dy- 
namics under gravity. 
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A. Spherical model: solutions in various approxi- 
mations 

Let us start from an initial profile 8i(x,ti) given at a red- 
shift Zi. We compute here various indicators, such as the 
density contrast, the divergence of the velocity field and 
the velocity field, at the time t, corresponding to present 
time so to z = 0. We assume that z 8 - ^ 1 so that transient 
terms have disappeared. Then, one can take into account 
only the faster increasing term. We also recall here the 
exact solution (Peebles, 1980). 

By integrating the total amount of matter inside a 
sphere (see Martel & Freudling, 1991), mass conservation 
is written 



q = Xi [l + S' i } 1 ^, 
with 



6'i 



y 2 6i(y)dy. 



(Al) 



(A2) 



Thus, the first order displacement M^ 1 ) is 

¥ r )(q) = ^{1-[1 + ^] 1 / 3 } (A3) 

The Lagrangian density contrast is written 



(A4) 



with x = q + (fli/flfi,,-)^ 1 ) + (</ 2 /</?,i)* (2) + 0(e 3 ). The 
calculation of easily leads to 



(2) 



(1) 



In the Eulerian case, we naturally obtain 



¥ le \x) 



1 (.i 

-o i 



V|/(2e) = Nj/(le)l I 



(A5) 



(A6) 



From this comparison to simulation for a initial power 
spectrum with little power at small scales it seems that 
second order Lagrangian approximation is a real improve- 
ment to first order approximation (Zeldovich) for on many 
aspects. It can describe higher density contrast, even if it 
is deficient to describe voids. And its validity after shell- 
crossing is smaller than with Zeldovich approximation. It 
gives surprisingly good results for the third, forth and fifth 
normalized moments of density distribution. L2 gives good 
results for velocity field, at least for its modulus. 

A.l. Spherical top hat model 

A. 1.1. The density contrast at the center of the perturba- 
tion 

Here, after having recalled the exact solution (Peebles 
1980), we compute, for various perturbative models, the 
density contrast 6 at the center of the perturbation. 6 does 
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not depend (before shell-crossing) on the considered pro- 
file. Spatial derivatives will be done with respect to the 
comoving coordinate x. 

Exact solution 

Once the value of Si = 6(0, ti) is given, the density contrast 
6 at time t can be analytically computed whatever the 
value of the density parameter £1. We recall here the exact 
solution (see Peebles 1980). 
Let us define the quantities 



F(Q) 



( 1 1 = 1 

2 1/S1 - 1 = chry- 1' 

_3_ Gt£ 

10 a ' 

1 1 _ 1 

2 1 - 1/S1 = 1 - cos j] ' 



n < i, 
si = i, 
n > i, 



6 C = -(1/fi,- - 1), 



si ^ i 
si = 1. 



(A7) 

(A8) 
(A9) 



Sli is the value of the density parameter at t = t;. It is 
given, as a function of SI, by the following expression: 



SI,; 



1 + Zi 



■SI. 



1 + Slz; 

With these definitions, 6(0) is 



(A10) 



8e(0) = { 



(<5 C - Si)/6 C 



F(ft)(ch9- 1) 
(<5 C - Si)/6 C 



F(tt)(l - cos( 



<5i < <5 C 
<5i > <$c 



(All) 



where 6 if the "proper conformal time" of the perturba- 
tion. It is determined by the implicit equation 



s h9-9 = \(6 c -6i)/6c\ 3 ' 2 G(n), 
e-sinO = \(6 c -6i)/6 c \ 3 ' 2 G(n), 



Si < 6 C 
Si > 6 C 



(A12) 



with 



G(S1) 



sh r) — T), 



5 a 
3 a,; 



3/2 



rj — sin T), 



SI < 1, 
SI = 1, 
SI > 1. 



(A13) 



Eulerian linear theory 

Equation (68), taken at first order in e, reads 

Mo) = — Si- 

9i, i 



(A14) 



Eulerian second order perturbation theory 
Equation (68) gives 



M0) = — 6i + 

9l,i 



9i 



9i, 



1 92_ 



6 



(A15) 



Therefore, the minimal density contrast that can be de- 
scribed in this approximation is 



c 3 / g 2 

Omin,E2 — 7 \ ^ n 

4 V g( 



Zeldovich approximation 



(A16) 



Mass conservation is here written as follows 



9i 



M0)= i + — p 

\ 9l,i 

with 

p = (i + Si)- 1/3 -i. 



Lagrangian second order perturbation theory 
Mass conservation reads 



(A17) 



(A18) 



M0)= (i + f-p + $-p> 

\ 9i, i 9ia , 



1. 



(A19) 



The smallest density contrast that can be reached in this 
approximation then is 



id 

4 92 



1. 



(A20) 



A. 1.2. The relationship -V.v/(aHS) = f(S) 



We assume here that Si(x,ti) is a continuously differen- 
tiable function for any x > 0. In the vicinity of the per- 
turbation center, the initial density contrast can then be 
expanded as Si(x,ti) = Si(0) + 0(x). A sphere of very 
small radius x is thus similar to a piece of homogeneous 
universe of density ~p\l + Si(0)]. So the proper velocity of 
an element of matter can be written: 



u(x) = a p H p x + 0(x 2 ), 



(A21) 



where a p and H p are respectively the expansion factor 
and the Hubble constant of this fictitious universe. The 
peculiar velocity then is 



v(x) = (a p Hp - aH)x + 0(x 2 ). 



(A22) 



The quantity we are interested in is written, in the limit 
x 



V.v 

aHS 



(0) 



3 ( a p Hp 



S \ aH 



(A23) 
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Exact solution 



and f2 is defined in Sect. 2. 



From Eq. (A23), the exact solution is simply (see also Zeldovich approximation 
Regos et al. 1989) 



V.v 

~~a~H6~ 



J E 



sh6»(sh( 



6 [Ht(ch9- l) 2 
3 
'8 



1 



Ht(l - cos( 



6 < 6 t 



6 > 6 t 



Equations (A17), (A18) and divergence of the first mem- 
ber of Eq. (70) give 



(A24) 



V.v 



iH6 



J Ll 



Jl + ^)- 1/3 ~l , 

6(1 + 6)- 1 ' 3 



(A31) 



with 



' 2 (chry - l) 3 
9 (sh r) — rj) 2 

0, 

2(1- COS T]) 3 

. 9 (r) — sin rj) 2 



ft < 1, 
ft = 1, 
ft > i. 



(A25) 



Lagrangian second order perturbation theory 
With Eqs. (A19) and (70), we easily obtain 

V.v' 



iH6 



J L2 



3 \ /l«L2 + £/2«£ 2 



6 J 1 + 6^2 + ^6? 



(A32) 



si 



L2 



where i^L2 is the solution (if it exists) of the following poly- 
The quantity Ht can also be written as a function of con- nom i a l 

formal time rj: 

' sh r)(sh rj — rj) 



Ht = t 



(ch rj — l) 2 ' 

2 
3' 

sin r)(r) — sin rj) 
. (1 — cos r)) 2 ' 



n < i, 
rt = i, 
ft > i. 



1 + «L2 + 

9i 



a +6) 



-1/3 



(A33) 



(A26) 



To obtain the quantity — V.v/(a_ffi5) as a function of 6, 
we need to suppress Si in the implicit equation (A12). We 
get 



(shfl - ff) 2 
(ch 6 - If 



2 1 + 6 

9 1 + ^turn 

2 1 + 6 



(l-cos6») 3 9 1 + (5 tu 



Eulerian linear theory 



6 < (Wn, 
^ > ^turn- 



(A27) 



^4.i?. Velocity field and density contrast for a given initial 
profile 

Here, we try to evaluate the density contrast and the veloc- 
ity field, as functions of radius, for a given initial profile 
6i(x). We recall the exact solution and give the various 
perturbative solutions. 

A.2.1. The density contrast 

Exact solution 

With the notations we use in this appendix, the exact 
solution can be written, for ft ^ 1, for respectively 6'i < 6 C 
and 6'; > 6 C 



Taking the divergence of expression (71) at first order in 
e reads, using Eq. (A14), 



V.v 

' aH6 



fi 



\l-6'i/6 c \ 3 [F(n)(ch6-l)y 



1 + 3- 



. 6i - 6'i 



(A28) Hr) = < 



3 sh 0(sh 9-9) 



2 (ch6»-l) 2 



El 



where fi is defined in Sect. 2. 

Eulerian second order perturbation theory 

Divergence of expression (71) gives, using Eq. (A15), 



(A34) 



\l-6'i/6 c \ 3 [F(Cl)(l-co S t 



Si - 6' 
For ft = 1, we have 



1 



3 sin ( 



2 (l-cos( 



V.v 

' aH6 



fi 



0E2 



E2 



92 f f 
—>52 - Jl 

91 



I2 
36 



(A29) 



where 6-^2 is the solution (if it exists) of the following 6(r) = < 
polynomial 



»E2 



v3 3g 2 



E2 



(A30) 



(sh6»-6») 2 (ch6»-l)- 



1 + 3(1 - Si /6'i) 



3sh6»(sh( 



2 (ch6»-l) 2 
sin6») 2 (l - cost?)" 3 



1 + 3(1 - 6i/6'i) 



3 sin 



2 (l-cos6») 2 
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for 6'i < 6 C and 6'i > S c respectively. This calculation is 
made in Lagrangian coordinates. In other words, the ra- 
dius r(t) = ax(t) where 6 is evaluated follows the motion. 
It is related to the initial radius r; = a; a;; by 



r = < 



3 ch 6 - 1 
10 ri 6 C - 6'i ' 



3 1 — cos ( 
:r;- 



10 6'i - 6 C 



6'i < 6 C 



6'i > 6 C 



(A35) 



The proper conformal time 6 is determined by the implicit 
equation (A12), in which Si has to be replaced by 6'i. 

Eulerian linear theory 

The density contrast is simply written 



<$ei(«) = — ^, 

9l,i 



Eulerian second order perturbation theory 
Equation (68) reads 

<5e2(«) = -^-Si 
9i,i 

9i,i J V 3 3 3 m 

92 /2 , 1 ,2 



(A36) 



(A37) 



Zeldovich approximation 

With the notations given in the beginning of this ap- 
pendix, we simply have 



M?)= i + 



9i 2 ( 1 9i d¥^\ 1 



1 



9i,i 1 ) \ 9i,i aq J 

This quantity is evaluated at x = q + /<?i,i) 1 ^ r< ~ 1 ' > - 
Lagrangian second order perturbation theory 
One easily has 



-l.(A38) 



M<?) = 1 + 



9 



1 §(1) g 2 §(2)~ 



0U rf 9 91 i aq / 



This quantity is evaluated at a; = q + (gi/gi i)^''- 1 - 1 + 

(Wff? 8 )* (2) - 



A. 2. 2. The velocity field 

We compute here the normalized velocity v/(Hri), where 
H is the Hubble constant and r\ is the typical radius of 
the considered fluctuation (so it exactly follows motion). 

Exact solution 

With the notations we use in this appendix, the exact 
solution can be written 



sh6>(sh( 



v(r) = I n Ht(di9 - l) 2 ' 
Hr\ ) r sin 6(6 — sin ( 



r\ Ht(l — cos ( 



6'i < 6 C 
6'i > 6 C 



(A39) 



This quantity is calculated at the radius r given by 
Eqs. (A35). 6 is still determined by the implicit equation 
(A12), but by replacing Si by 6'i. 

Eulerian linear theory 

At first order, equation (71) reads 



v(x) 



1 r g 



with r = ax. 



3 Hri gi 



1 x' 
— 6 ,;, 



(A40) 



Eulerian second order perturbation theory 
At second order, equation (71) reads 



v(x) 



Er x 



1 r 



J E2 



3 Hn 
fflffl 

" a 2 
9ia 



9i 1 92 w2 

9i,i 3 gf ti 

cl c 2 „,2 

i°i - i 



(A41) 



with r = ax. 

Zeldovich approximation 

Taking Eq. (70) at first order simply leads to 



v(q) 



Er x 



= _EL_£L*(i) j 

li Hr x gi ti 



(A42) 



This quantity is evaluated at x = q + (gi/ gi^)^ 1 ^ ■ 
Lagrangian second order perturbation theory 
At second order, Eq. (70) reads 

v(q) 



Er x 



= JL- 1L§(1) + 1^§(2) 

L2 #>i \g h i g\ i 



(A43) 



This quantity is evaluated at x = q + (gi/gi i)^ 1 " 1 ^ + 

(Wff?,.0* (2) - 
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